function [value,policy,probs] = solve_bellman(S,T,lambda,phi,y_high,y_low,p)

% optimal number of kids 
n_star = length(S)-2; % S starts with n=0 and ends with n=n_star+1

% fertility path
eta = (1-lambda) * phi;

% initialise value function
value = zeros(3,n_star+1,T);
policy = zeros(3,n_star+1,T);
value(:,:,T) = ones(3,1) * S(1:end-1);

% loop backwards over time periods (we know values for T)
for t = T-1:-1:1
    
    % ZERO CURRENT INCOME
    
    % loop over current n
    for n=0:n_star 
        % corresponding index of value array (n=0 is element 1)
        ix = n+1; 
        % if target has been reached, do not get pregnant
        if n == n_star 
            value(1,ix,t) = S(ix);
            policy(1,ix,t) = 0;
        % if target has been reached, get pregnant 
        else
            % success with probability eta(t)
            value(1,ix,t) = eta(t) * value(1,ix+1,t+1) + (1-eta(t)) * value(1,ix,t+1);
            policy(1,ix,t) = 1;
        end        
    end
    
    % LOW INCOME
    if n_star > 0 
        v_preg = eta(t) * value(1,2,t+1) + (1-eta(t)) * value(1,1,t+1);
    else 
        v_preg = 0;
    end
    v_work = y_low + p * value(3,1,t+1) + (1-p) * value(2,1,t+1);
    value(2,1,t) = max(v_preg,v_work); % bellman
    policy(2,1,t) = (v_preg > v_work);
    
    % HIGH INCOME
    if n_star > 0
        v_preg = eta(t) * value(1,2,t+1) + (1-eta(t)) * value(1,1,t+1);
    else
        v_preg = 0;
    end
    v_work = y_high + value(3,1,t+1);
    value(3,1,t) = max(v_preg,v_work); % bellman
    policy(3,1,t) = (v_preg > v_work);
end


end
